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We present a generic approach for treating the effect of nuclear motion in the high-order harmonic 
generation from polyatomic molecules. Our procedure relies on a separation of nuclear and electron 
dynamics where we account for the electronic part using the Lewenstein model and nuclear motion 
enters as a nuclear correlation function. We express the nuclear correlation function in terms of 
Franck-Condon factors which allows us to decompose nuclear motion into modes and identify the 
modes that are dominant in the high-order harmonic generation process. We show results for the 
isotopes CH4 and CD4 and thereby provide direct theoretical support for a recent experiment [Baker 
et al., Science 312, 424 (2006)] that uses high-order harmonic generation to probe the ultra- fast 
structural nuclear rearrangement of ionized methane. 
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I. INTRODUCTION 

In the past decade several experiments have demon- 
strated the potential for high-order harmonic generation 
(HHG) to probe molecular structure. These experiments 
include studies of orbital structure [TH3] , nuclear dynam- 
ics [HE] and more recently coupled electronic and nuclear 
dynamics [7J [8] . The interplay between nuclear motion 
and HHG has been subject to several theoretical stud- 
ies in the case of H^ and H2 |9HT2] and other diatomic 
molecules [13] where nuclear motion is simple in the sense 
that there exists only one mode of vibration (along the 
internuclear axis). So far, only limited theoretical work 
deals with the influence of nuclear motion on HHG from 
non-linear molecules [141 US] • 

Five years it was predicted that HHG may be used to 
probe the fast nuclear motion in H2 by comparing the 
HHG spectrum to that of the isotope D2 [10 , and this 
was soon after confirmed in a pioneering experiment [5]. 
The D2 /H2 ratio of the harmonic spectra is an increasing 
function of the harmonic order which can be understood 
from the basic idea that the harmonic order is associated 
with the time the electron spends in the continuum from 
initial ionization to recombination, r. The larger r, the 
higher the harmonic order. At 800 nm and intensities 
around 10 14 W/cm 2 typical r's are in the 1-2 fs regime 
and the reason for the observed increase is the faster nu- 
clear motion in the lighter isotope H2 , leading to a smaller 
overlap of the nuclear wave packets in the recombination 
step. In the experiment [5] similar results were reported 
for HHG from methane isotopes (CD4/CH4). 

The current work addresses the effect of nuclear motion 
on HHG for arbitrary, linear or non-linear, polyatomic 
molecules. As we will see, the effect is generally incor- 
porated through a nuclear correlation function describing 
the nuclear wave packet dynamics in the molecular cation 
from initial ionization to recombination. This function is 
highly demanding to determine. To solve the problem, we 
first express the nuclear correlation function in terms of 



Franck-Condon (FC) factors (denned as the square of the 
overlap integral between the vibrational wavefunctions of 
the neutral molecule and the molecular cation, in their 
respective electronic states) and the accompanying time- 
dependent phases caused by the vibrational excitations. 
Second, we consider these in the harmonic approxima- 
tion, where they may be calculated using standard ap- 
proximations and technology from quantum chemistry. 
The theory is exemplified by calculations on CH4 and 
CD4. Methane has nine modes for vibrational relaxation, 
but the model used here allows us to identify the two 
most important vibrational modes which turn out to pro- 
mote aTdH C2V nuclear geometry reconfiguration [see, 
e.g., Ref. [16 for a discussion of point group symmetry]. 
We thereby prove a conjecture put forward in [5]. 

This paper is organized as follows. In Sec. [Hj we 
present the theory. The electronic part is well-known 
from other works, so we focus on the analysis of the vi- 



brational motion. In Sec. Ill, we present results on CH4 
and CD4. Section [IV| concludes. In Appendix |A| we dis- 
cuss the calculations of the FC factors. Atomic units 
[h = e = rn e = clq = 1] are used throughout this paper. 



II. THEORETICAL MODEL AND 
COMPUTATIONAL DETAILS 

We wish to provide a simple model that in a general 
and transparent way isolates the role of nuclear motion in 
the typical experiment on HHG from molecules in the gas 
phase. We treat electronic dynamics in the Lewenstein 
model [17] following the implementation of Ref. [4 , but 
adapt a few improvements to the latter as detailed below. 

The observable quantity is the HHG spectrum, which 
we calculate from the Fourier transform of the dipole 
velocity (v(t)) [T8] , .i.e., 
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where the laser field driving the process is non-zero only 
in the time interval [0,T]. Since the early 1990's it has 
been an ongoing discussion whether the HHG spectrum 
should be calculated from the (pulse limited) Fourier 
transform of the dipole moment, the dipole velocity or 
the dipole acceleration. For relatively long and weak 
laser pulses the final result is independent of the choice: 
Up to a well-known frequency dependent factor one can 
interchange the dipole velocity with the dipole moment 
or the dipole acceleration, since the appropriate bound- 
ary terms vanish almost exactly. For short and intense 
laser pulses this is no longer true, and the use of a wrong 
form leads to an unphysical background in the HHG spec- 
trum [TO] [20] . We have tested our calculations using the 
different forms, and we find good agreement between the 
velocity and the acceleration form, whereas the length 
form differs considerably. These findings are consistent 
with previous studies [21] , and so we restrict our analysis 
in this paper to the velocity form results. 

To evaluate the dipole velocity (v(t)) = 
in Eq. we apply the Born-Oppenheimer approxima- 
tion, the strong-field approximation and single- active- 
electron model. We further freeze all orbit als except 
the highest occupied molecular orbital (HOMO) and 
consider the orientation of the nuclei as fixed during 
the short high-harmonic generating femtosecond pulse 
F(t) = F(t)e of linear polarization e. Then introduc- 
ing the molecular orbital of the active electron evalu- 
ated at the nuclear equilibrium configuration, ^$(t) = 
ifjo(r) ex.p(il p t) (I p is the adiabatic ionization potential 
of the molecule), the ground vibrational state of the neu- 
tral molecule, Xi,o> and using a product of a Volkov wave, 
ifrk , and the uth vibrational state of the molecular ion, 
Xf,u fc> r the propagator we have 

|*(t)> = |*„Xi,o(t)> 

dt' f d s kY J \^Pl{t)x s At)) 

x {i>l(t') X fAt')\F(t')e • r\Mt')XiAt')) (2) 
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is the vibrational autocorrelation function (Fig. [T]) with 
l(Xi>IXi,o)| 2 FC factors (see Appendix |a| , e v the vibra- 



tional energy, 

v rec (k) = /c(2tt)" 3/2 J d 3 rexp[-ik • r]i/j (r) 

= fe0o(fe), (5) 
dion(k) = e • (27r)" 3/2 J d 3 rex.p[-ik • r}ri/j (r) 

= ie>V k <l>o(k), (6) 
5(fc, t, if) = jT dt"[(k + A(t")f/2 + J p ], (7) 

and A(t) the vector potential of the laser field. 

The electronic part of Eq. (|3| has appeared many times 
since the seminal paper [17] and clearly points out the 
three essential steps of HHG process: The electron ion- 
izes to the continuum at time t' with probability ampli- 
tude F{t f )d[ on (k + A (£')). It then propagates in the field 
until time t acquiring a phase factor S(fe, t, t') and recom- 
bines with a probability amplitude v* ec (k + A(t)). Due 
to vibration this product of amplitudes is weighted by a 
nuclear factor C(t — t r ) [TO] . 

To evaluate the electronic part of Eq. ([3| we write the 
molecular orbital of the active electron in the molecular 
fixed (MF) frame as a linear combination of Gaussian 
orbit als that we find using the GAMES S quantum chem- 
istry code [22 . Our basis choice is TZV with polarization 
functions. The expansion of the molecular orbit als in 
terms centered on the atoms allows us to evaluate the k- 
integrals of Eq. ([3| within the improved stationary phase 
method [23-25 , where we include electron trajectories 
leading from one atomic center to another. Dipole ve- 
locities for different orientations of the molecule can be 
obtained by applying the Euler rotation operator [4j [26] 
to the molecular fixed wave function. 

III. RESULTS AND DISCUSSION 

We now turn to an application of the theory described 
above. We consider HHG from CH4 and CD4. To cal- 
culate the harmonic yields, we first determine the elec- 
tronic structure and the FC factors. Using GAMESS we 
find the following configuration of the molecular orbitals 
(lai) 2 (2ai) 2 (3t 2 ) 6 and thus we have six electrons dis- 
tributed among three degenerate HOMOs (see Fig. [2| 
that could all contribute appreciably to the harmonic 
yield. The adiabatic ionization potential (I p — 12.92 
eV) is estimated by comparing the total energy of the 
methane to that of the relaxed methane ion. The effec- 
tive ionization potential is then obtained by adding the 
additional shifts (e^) due to vibrational excitation. Our 
FC analysis shows that only two normal modes are ex- 
cited when CH4 (CD4) ionizes. These are an E symmetry 
mode (u ^1295 cm" 1 for CH+; uo ^920 cm" 1 for CD+) 
that brings the molecule towards a plane and an Ai sym- 
metry mode (a; ~2766 cm -1 for CHj; uj ~1960 cm -1 for 
CD^) that correspond to changes of the C— H (C— D) 
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FIG. 1: (Color online) Simple sketch of the effect of vibra- 
tion on HHG. When the molecule ionizes the laser launches 
a vibrational wave packet that evolves on the ionic Born- 
Oppenheimer surface. The overlap of this wave packet and 
the initial vibrational state [see Eq. Q] weights the dipole 
velocity [see Eq. (p}] and hence the HHG signal. 
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bond lengths. In Fig. [3] we show the one-dimensional FC 
factors for both modes (calculated from integrals over 
a single mode as shown in Appendix [A| along with in- 
sets indicating the nuclear rearrangements related to each 
normal mode. These modes drive the molecular ion into 
the relaxed C2 V symmetry. 

With the molecular structure at hand, i.e., molecular 
orbitals and FC factors, we should, in principle, calcu- 
late the molecular dipole velocity as given in Eq. ([3) for 
different orientations of the molecule and for each of the 
degenerate orbitals and average the resulting dipole ve- 
locities [see Ref. [27]]. However, due to the high symme- 
try of the methane molecule there is only a small angular 
dependence of the harmonic yield and further the yield 
is dominated by the contribution from a single orbital 
that couples strongly to the linearly polarized laser field. 
We have verified this by calculations for various orien- 
tations not shown here. Consequently, the averaging is 
redundant and the results presented in this paper are car- 
ried out for one fixed orientation of the molecule and in- 
cludes only the yield of the dominant orbital (cf. Fig. [2]). 
Another technical detail is that we limit the t' integral 
in Eq. ([3| to times when r = t — t' is smaller than 0.65 
times an optical cycle. We do this to have a one-to-one 
correspondence between the time the electron spends in 
the continuum (from the instant ionization, t f , to the mo- 
ment of recombination, t) and the harmonic order, i.e., 
the energy released when the electron recombines [TP] . 

Figure [4] shows the harmonic spectra of CH4 and 
CD4. We use a trapezoidal shape for the vector poten- 
tial A(t) = A(t)e x , with two optical cycles linear turn-on 
and turn-off and three cycles of constant amplitude cor- 
responding to peak intensity of 2 x 10 14 W/cm 2 . The car- 
rier wavelength is 775 nm. The curves shown on the fig- 



FIG. 2: (Color online) Isocontour plots of the degenerate HO- 
MOs of methane. We calculate HHG yields with the linear 
laser field polarization along the x axis and include only the 
contribution from HOMOl. 
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FIG. 3: (Color online) Franck-Condon factors for the E mode 
(panel (a)) and Ai mode (panel (b)) of CH4 (dashed) and 
CD J (full). These two dominating modes will drive the 
molecule from the T^ symmetry into the relaxed C2V sym- 
metry upon ionization as conjectured in [5]. 
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FIG. 4: (Color online) Harmonic spectra for CH4 and CD4 
using a 775 nm linearly polarized laser of intensity 2 x 10 14 
W/cm 2 and the trapezoidal envelope detailed in the main 
text. 

lire exhibit the typical characteristics of harmonic spec- 
tra, viz., an exponential drop-off at low harmonic orders 
followed by a plateau with a cutoff around the harmonic 
order 33 in agreement with the cutoff formula [T7] . 

It is hard to see the difference in the spectra for the 
two isotopes, but if we integrate the spectra in an interval 
around each (odd) harmonic order and plot the ratio of 
these numbers for the isotopes, we get the result shown 
in Fig. [5] by the dash-dotted curve. For comparison, we 
also plot the CD4/CH4 ratio of the nuclear correlation 
functions (see Eq. (Ek) if we include only the dominant 
E symmetry mode (Ce), both the E and the Ai symme- 
try mode (Ce+aJ an d if we further include two slightly 
excited T2 symmetry modes (Ce+Ai+t 2 )- We see that 
the E+Ai motion accounts for the major effect of nu- 
clear motion. Further, the nuclear correlation functions 
are increasing monotonically and hence in our model the 
oscillatory structure of the HHG ratio arises from the 
combined electron- nuclear dynamics. 

The ratio predicted by the current theory underesti- 
mates the slope of the ratio as compared to measure- 
ments [5] . To understand this deviation we refer to more 
detailed calculations carried out on the simpler systems 
H2 and D2 . For these systems we may clearly understand 
the consequences of the approximations made here, and 
as such attribute the disagreement to two factors. First, 
we expect that the harmonic approximation, used here to 
retrieve the FC factors, yields nuclear factors, C(t — 
with a too low ratio between the isotopes since stretching 
of the molecule is underestimated when the asymmetry 
of the potential is not taken into account. To substanti- 
ate this conjecture we have checked the case of H2 and 
D2, where we can compare the nuclear correlation func- 
tions resulting from FC factors based on the harmonic 
approximation and a more accurate Morse potential [28 , 
respectively and our reasoning is validated (see Fig. |6t. 
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FIG. 5: (Color online) Curves displaying the ratio of the 
harmonic spectra from Fig. [4] (dash-dotted) and of the 
CD4/CH4 nuclear correlation functions (Eq. Q) including 
normal modes of symmetry E (dashed), E and Ai (full) and 
E and Ai and T 2 (dotted). 
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FIG. 6: (Color online) Ratio of the D2/H2 nuclear correlation 
function (Eq. Q) for the harmonic approximation (solid) and 
the Morse potential (dashed). 



Second, we have not included any coupling of nuclear 
and electron dynamics. In HHG from the isotopes H2 
and D2 such coupling is known to result in a dynamic 
two-center interference effect that leads to a higher ratio 
in the D2/H2 spectrum 0. 

The nuclear correlation function is expressed in terms 
of the FC factors (Eq. Q). To this end we note that 
for systems where the weight on the different FC factors 
are effectively experimentally adjustable by laser pulse 
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preparation one may control the nuclear correlation func- 
tion. This is for instance the case in H2 [29] and we have 
checked that by populating only FC factors (0, v — 0) 
and (0, v = 18) the nuclear correlation function oscillates 
with a period of 1.56 fs. Consequently, the nuclear cor- 
relation function can be modulated within less than half 
an optical period of the typical HHG driving field (1.33 fs 
for an 800 nm field) and thereby drive enhancement and 
suppression of certain electron trajectories and associated 
harmonic orders. Control of the relative strength of the 
harmonics is for instance useful for attosecond pulse gen- 
eration [30] and the control scheme briefly discussed here 
is based on the intrinsic structure of the molecule. Fi- 
nally, since the nuclear correlation function is expressed 
in terms of one-dimensional FC factors, we can identify 
the important part of the nuclear dynamics by study- 
ing the changes of the nuclear correlation function as we 
gradually add the different normal modes and by includ- 
ing the minimal amount of modes in the full HHG calcu- 
lation we save CPU time. 



Appendix A: Computation of the nuclear factor of 
Eq. 

We are interested in evaluating Franck- Condon inte- 
grals, i.e., (xf,v\Xi,o), between a vibrationally-cold initial 
state (Xi,o) and a vibrationally-excited final state (xj>) 
where the subscript v denotes the excitation level. 

In the normal-coordinate representation within the 
harmonic approximation, the vibrational wavefunctions 
describing the initial state (i.e., the neutral molecule) and 
the final states (i.e., the ion), respectively, are expressed 
as 



1 



XiM) = (detr'/Tr^)* exp ( —Q'Vq' ), (Al) 



and 



X/,„(Q) = (det 3 exp ( -^TQ 



^(2^n 3 \)- 1/2 H n3 (T^ 2 Q), (A2) 



IV. CONCLUSION AND OUTLOOK 

In conclusion, we have followed [10] and applied the 
Lewenstein model [17 to molecules with moving nuclei. 
We assume that the electronic and nuclear part separate. 
The electronic part is then treated conventionally. For 
the nuclei, however, we relate the vibrational autocorre- 
lation function to FC factors [see Eq. Q] and associated 
dynamical phase factors. For some polyatomic molecules 
the FC factors and energies are available in the literature, 
and the vibrational part of the theory can be determined 
directly. 

In the cases where FC factors and vibrational ener- 
gies are not available, we discuss how to perform a nor- 
mal mode analysis and calculate the intensity of the FC 
transitions in the harmonic approximation. The model 
covers any polyatomic molecular system where the elec- 
tron dynamics is reasonably described within the single- 
active-electron picture. The theoretical and computa- 
tional models involved are fairly standard within strong- 
field physics and computational chemistry. Finally, with 
the present theory, we identify the vibrational modes in- 
volved in the ultrafast rearrangement of the nuclei in the 
CH4/CD4 system, and obtain qualitatively agreement 
with the measurements [5]. 
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where Q = (Qi,Q2 5 *"' iQN m ) is a normal coordinate 
vector, r is a diagonal matrix with elements Fjj = h/ujj 
where Uj is the vibrational frequency of mode j, N m is 
the number of vibrational degrees of freedom, the index 
n = (ni,n2, • • • , ^iv m ) is a vector of vibrational excita- 
tions (n is related to v through the function A that orders 
the excited modes according to ascending energy, i.e., 
v{n) = A(n)) and H nj is the n^th Hermite polynomial. 
For a vibrationally-cold state, n = (0i, O2, • • • , 0jv m ) = 0. 
The normal coordinates of the initial state (Q 1 ) and 
the final state (Q) are related by a simple transforma- 
tion [51] , i.e., 



Q' = JQ + A, 



(A3) 



where J is the so-called Duschinsky matrix and the vector 
A expresses the geometry change in the final state. The 
J matrix reflects the mapping of the normal coordinates 
of the initial-state onto those of the final-state. 

The multi-dimensional Franck-Condon integral, within 
the harmonic approximation, reads 

(Xf,n I Xi,o) = N J dQi . . .dQ Nrn H ni (T 1 Q 1 ) . . . 



x^ m (r 



x exp 
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where the normalization factor is given by 

Nrr 
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When evaluating the integrals, a considerable simpli- 
fication is introduced b y as suming that the off-diagonal 
elements in J (see Eq. (|A3[)) are very small, i.e., 



Jj,iQi + Jj^Q2 



Jj,N m QN n 



A, 



■A,. 



(A6) 



Accordingly, the multi-dimensional Franck-Condon in- 
tegral reduces to a product of one-dimensional integrals, 
i.e, 



The Ansbacher recurrence relations are used to ob- 
tain one-dimensional FC integrals [32 . Computing FC 
integrals involves calculation of the equilibrium geome- 
tries and normal modes for the neutral and the ionized 
molecules. These were obtained from calculations using 
the hybrid density functional B3LYP level of theory in 
conjunction with the triple-^ valence basis set as imple- 
mented in Gaussian 33 . 



(Xf,n I Xi,o) = [I / dQjH^iTjQj) 



x exp 



(A7) 
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